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We employ the full FENE-P model of hydrodynamics of a dilute polymer solutions to 
derive a theoretical approach to drag reduction in wall bounded turbulence. We recapture 
the results of a recent simplified theory which derived the universal Maximum Drag 
Reduction (MDR) asymptote, and complement that theory with a discussion of the 
cross-over from the MDR to the Newtonian plug when the drag reduction saturates. 
The FENE-P model gives rise to a rather complex theory due to the interaction of 
the velocity field with the polymeric conformation tensor, making analytic estimates 
quite taxing. To overcome this we develop the theory in a computer-assisted manner, 
checking at each point the analytic estimates by Direct Numerical Simulations (DNS) of 
viscoelastic turbulence in a channel. 



1. Introduction 

The onset of turbulence in fluid flows is accompanied by a significant increase in the 
drag Lumley (1969), Sreenivasan & White (2000). This drag poses a real technological 
hindrance to the transport of fluids and to navigation of ships. It is interesting therefore 
that the addition of long chain polymers to wall-bounded turbulent flows can result in a 
significant reduction in the drag. The basic experimental knowledge of the phenomenon 
had been reviewed and systematized by Virk (1975); the amount of drag depends on 
the characteristics of the polymer and its concentration, but cannot exceed a universal 
asymptote known as the "Maximum Drag Reduction" (MDR) curve which is independent 
of the polymer's concentration or its characteristics. When the concentration is not large 
enough, the mean velocity profile as a function of the distance from the wall follows the 
MDR for a while and then crosses back to a Newtonian-like profile, cf. Fig. 1 left. 

Recently the nature of the MDR and the mechanism leading to its establishment were 
rationalized, using a phenomenological theory in which the role of the polymer confor- 
mation tensor was modeled by an effective viscosity L'vov et. al. (2003). It turned out 
that the effective viscosity attains a self-consistent profile, increasing linearly with the 
distance from the wall. With this profile the reduction in the momentum flux from the 
bulk to the wall overwhelms the increased dissipation that results from the increased vis- 
cosity. Thus the mean momentum increases in the bulk, and this is how drag reduction 
is realized. In De Angelis et. al. (2004) it was demonstrated by DNS that Navier-Stokes 
flows with viscosity profiles that vary linearly with the distance from the wall indeed 
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Figure 1. Left panel:Mean normalized velocity profiles as a function of the normalized distance 
from the wall during drag reduction. The data points from numerical simulations (green circles) 
and the experimental points (open circles) represent the Newtonian results. The red data points 
(squares) represent the Maximum Drag Reduction (MDR) asymptote. The dashed red curve 
represents the theory discussed in the paper which agrees with the universal MDR. The arrow 
marks the crossover from the viscous layer to the Newtonian log-law of the wall. The blue filled 
triangles and green open triangles represent the crossover, for intermediate concentrations of the 
polymer, from the MDR asymptote to the Newtonian plug. Right panel: Mean velocity profiles 
for the Newtonian and for the viscoelastic simulations with Re T — 125. Solid line: Newtonian. 
Dashed line: Viscoelastic. The straight lines represent the classical von-Karman log-law. The 
coordinate y + is the distance from the wall in Pradtl'units y + = v*L/v, where v't = p L. 

show drag reduction in close correspondence with the phenomenon seen in full viscoelas- 
tic simulations. The aim of this paper is to complement the simplified theory with a 
derivation of the same results on the basis of a full viscoelastic model of the hydrody- 
namics of dilute polymer solutions. Such a derivation forsakes some of the simplicity of 
the phenomenological theory, but on the other hands it clarifies the role of the polymer 
conformation tensor in furnishing viscosity- like contributions. The relative complexity of 
the statistical theory is handled by employing suitable approximations for the leading 
terms and by computer assisted estimates of competing terms, allowing us at the end to 
provide a compact theory with sharp predictions. In addition to lending further support 
to the model of L'vov et. al. (2003), we will offer a discussion of the non-universal satu- 
ration of drag reduction as a function of concentration, length of polymer and relaxation 
time in various conditions of experimental interests. 

In Sect. 2 we consider the FENE-P model of viscoelastic flows and review shortly 
the evidence for drag reduction in DNS of this model. In Sect. 3 we employ the FENE-P 
model to derive a statistical theory of drag reduction in wall bounded turbulence. In Sect. 
4 we use the theory to predict the cross-over from the universal MDR to the Newtonian 
plug when the conditions differ from those necessary for attaining the MDR. In Sect. 5 
we present a summary and a discussion. 

2. Equations of motion for viscoelastic flows and drag reduction 

Viscoelastic flows are represented well by hydrodynamic equations in which the effect 
of the polymer enters in the form of a "conformation tensor" Rij (r, t) which stems from 
the ensemble average of the dyadic product of the end-to-end distance of the polymer 
chains Bird et. al. (1987), Beris & Edwards (1994). A successful model that had been 
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employed frequently in numerical simulations of turbulent channel flows is the FENE-P 
model Bird et. al. (1987). Flexibility and finite extensibility of the polymer chains are 
reflected by the relaxation time r and the Peterlin function P(r, t) which appear in the 
equation of motion for Rij : 

dR a p _ du a du 

-\[P{r,t)R a p - pl8 aP ] , (2.1) 

P(r,t) =(p 2 m -p 2 )/(pl-R J ,) • (2.2) 

In these equations and p\ refer to the maximal and the equilibrium values of the 
trace i? 77 . Since in most applications p m S> po the Peterlin function can also be written 
approximately as P(r,t) w (1/(1 — ai? 77 ) where a = p~f. In its turn the conformation 
tensor appears in the equations for fluid velocity u a (r, t) as an additional stress tensor: 

"" a + (u 7 V 7 )u Q = -V a p + v s V 2 u a + V 7 T Q7 + F a , (2.3) 



dt 

-P(r,t) 



T a0 (r,t) = ^ 



2 Rap{r,t)-8 a p 
Po 



(2.4) 



laf} ^ JLl_ Ral3 , (2.5) 



Here v s is the viscosity of the neat fluid, F a is the forcing and v p is a viscosity parameter 
which is related to the concentration of the polymer, i.e. v p jv s <~ c p where c p is the volume 
fraction of the polymer. We note however that the tensor field can be rescaled to get rid 
of the parameter a in the Peterlin function, R a p = aR a /3 with the only consequence 
of rescaling the parameter v p accordingly. Thus the actual value of the concentration is 
open to calibration against the experimental data. Also, in most numerical simulations, 
the term Pp^ 2 R is much larger than unity and the unity tensors in (2.1) and (2.4). 
Therefore, in the theoretical development below we shall use the approximation 

VpP_ 

T pI' 

These equations were simulated on the computer in a channel or pipe geometry, reproduc- 
ing the phenomenon of drag reduction in experiments. The most basic characteristic of 
the phenomenon is the increase of fluid throughput in the channel for the same pressure 
head, compared to the Newtonian flow. This phenomenon is demonstrated in Fig. 1 right 
taken from De Angelis et. al. (2003). As one can see the simulation is limited compared 
to experiments; the Reynolds number is relatively low, and the MDR is not attained. 
Nevertheless the phenomenon is there, and we will be able to use the simulation to asses 
and support the theoretical steps taken in the theoretical development. We should note 
that in the simulations one needs to add a small artificial viscosity to Eq. (2.1), in the 
form of a term n/S.R a p. This is done solely for taming numerical instabilities, and has 
very little consequence on the theory presented below, where it is not taken into account. 



3. The derivation of the MDR 

3.1. Exact balance equations 
Consider the fluid velocity U a (r) as a sum of its (time) average and a fluctuating part: 

U a (r,t) = V a (y)+u a (r,t) ,V a (y) = (U a (r,t)) . (3.1) 

For a channel of large length and width all the averages, and in particular V a (y) — ► 
V(y)S a y, are functions of y only. The objects that enter the theory are the mean shear 
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S(y), the Reynolds stress W(y), the kinetic energy K{y) and the mean conformation 
tensor (R(y)); the first three are defined respectively as 

S(y) = dV(y)/dy , W (y) = - (u x u v ) , K (y) = (\u\ 2 )/2 . 

Taking the long time average of Eq. (2.3), and integrating the resulting equation along 
the y coordinate produces and exact equation for the momentum balance: 

W + vS+^(PR xy ){y)=p'{L-y) . (3.2) 

T 

where p' is the pressure gradient and L is the mid-channel height. The RHS is simply 
the rate at which momentum is produced by the pressure head, and on the LHS we have 
the Reynolds stress which is the momentum flux, the viscous dissipation of momentum, 
and the rate at which momentum is transferred to the polymers. Near the wall it is 
permissible to neglect the term p'y on the RHS for Next we want to derive the 

balance equation for energy. To this aim we need to take into account the extra energy 
dissipation due to transfer of energy from the velocity field to the polymers. This energy 
dissipation, denoted as e p , in the limit of validity of (2.5) can be exactly computed to be 

e p = ^{P 2 TvR aP ) , (3.3) 

where TrR a p = R xx + R yy + R zz . As stated in Eq. (3.1), we split the velocity field into 
its mean and the fluctuation, and consider separately the balance equation for the mean 
energy V 2 and for the turbulent energy (u 2 ). The former yields an equation identical to 
(3.2) but multiplied by S: 

WS + vS 2 + ^- (PR xy ) (y)S = p'LS . (3.4) 

T 

On the other hand the balance equation for the turbulent energy cannot be written 
exactly. To understand how to write it, we need to identify first the meaning of the third 
term on the LHS of Eq. (3.4). To gain insight on how to represent properly this term we 
consider how to estimate correlations containing the tensor Rij. Using equation (2.1), we 
obtain on the average: 

d 1 

— (uyRij) = --(PRij) + (RikdkUj) + (RjkdkUi) . (3.5) 

3.2. Closure approximations 

To proceed, we employ a 1-point closure approximation, keeping only terms that sur- 
vive this closure. In this approximation, the only terms (dkUjRik) which survive are 
proportional to S, i.e. 

(RikdkUj) — aiS(R iy }5 jx (3.6) 
where ai are constant of order 1. For i = x,y, using (3.6), equations (3.5) give: 

-(P)(R XX ) =2a x (R xy )S , (3.7) 

T 

1 - { P)(R xy )=a y{ R yy )S. (3.8) 

To test the quality of this closure we check each term by DNS, comparing between the 
exact and approximated forms. A description of the DNS simulations is reported in De 
Angelis et. al. (2003). 

In Figs. (2) we show the comparison of the 1-point closure approximation, using nu- 
merical simulations with a x = 1 and a y — 0.45. The overall comparison is reasonably 
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Figure 2. Left panel: quality of the 1-point closure for equation (3.7) as obtained by direct 
numerical simulation for a turbulent channel flow, see De Angelis et. al. (2003). The continuous 
line refers to the term 2a x R xy S with a x = 1 while the dashed line to the term {P){R xx )/t. 
The the two terms are of the same order of magnitude within the buffer region. Right panel: 
comparison of the 1-point closure approximation for equation (3.8). The continuous line refers 
to the term a y R yy S with a y — 0.45, while the dashed line to the term (P){R xv )/t. The the two 
terms are of the same order of magnitude within the buffer region. 



good within the buffer layer, i.e. in the region where the polymers enhance the mean flow 
with respect to the Newtonian case. We can thus safely proceed and learn from Eq. (3.7) 
that in 1-point closure we can write the third term on the LHS of Eq. (3.4) as 



1 ±(PR xy )(y)S=^(Pf(R xx ) ; 



(3.9) 



as such, it is one of the three terms appearing in the trace in Eq. (3.3). As this term is 
already exhausted in the balance for the mean energy V 2 , this leaves for the turbulent 
energy balance equation the other two terms. We stress here that this conclusion rests 
entirely on the numerical value a x = 1. We propose that this is essentially exact, even 
though at this stage we present it on the basis of numerical comparisons. In future work 
we will derive this important fact from first principles. 

We can therefore write now the balance equation for the turbulent energy. The rate 
of production is simply W(y)S(y). This production of energy is balanced by the energy 
dissipation ^(|Vu| 2 ) and the two terms in the energy transferred to the polymer (3.3). 
The energy dissipation is estimated differently in the viscous layer near the wall and in 
the bulk of the turbulent flow, namely 

R 3/2 



K 

av— + b- 

y z 



A 



y 



+ ^(Py((Ryy) + (R ZZ )) = WS , 



(3.10) 



where the first term reflects the smoothness of the field in the viscous layer, and the 
second term is a standard Kolmogorov-type estimate of the energy flux as the typical 
energy at distance y from the wall over the typical eddie turn-over time yj \JK (y). These 
estimates are identical to the Newtonian case. The new (and somewhat difficult) term is 
the third term on the LHS, where A is a constant to be computed. 

To evaluate this last term, we refer again to our DNS results to assess the relative 
importance of (R yy ) and {R zz ): these objects are very close to each other throughout the 
region of concern in the channel. We therefore can keep only one of them by redefining 
consequently the constant A: 



2 K . 
a 2 u— + b 

y 



K 3/2 A 2 



y 



+ 



2-r(Pr(Rvv) 



ws 



(3.11) 
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3.3. Reminder: the von-Kdrmdn law 

At this point it is useful to remind to the reader the derivation of the von-Karman log law 
of the wall in the Newtonian case, see Pope (2000). The terms including the conformation 
tensor in Eqs. (3.2) and (3.11) are not present in this case. Outside the boundary layer 
the viscous terms proportional to v are neglected, and we therefore get from Eq. (3.2) 
the result that W is constant, W w p'L, and from Eq.(3.11) that bK 3 ^ 2 /y ~ WS. Now 
one introduces the experimental knowledge that K and W are proportional to each other 
outside the viscous boundary layer: 

W = c 2 N K . (3.12) 

This equation can be presented as a rigorous inequality W < K (i.e. c?f is at most unity), 
and it can be justified phenomenologically as a result of the scale invariance. We will take 
it, and its visco-elastic counterpart 

W = <? V K, (3.13) 

as given experimental inputs. Using (3.12) in the relations obtained for the Newto- 
nian case, we find immediately S oc y/pTL/y which integrates to the log-law of the wall 
V = logy/nx + B where the von-Karman constant kk and the intercept B are phe- 
nomenological constants. 

3.4. Derivation of the MDR 

We anticipate that in the viscoelastic case the dominant terms in the balance equations 
(3.2) and (3.11) will be the terms including the polymer conformation tensor. To see that 
this is indeed the case, we analyze the yy component of Eq. (3.5). It reads 

^(UyRyy) = ~{PRyy) + 2(R yk d k u y ) . (3.14) 

In 1-point closure this equation gives the absurd result {PR yy ) — which means that in 
this case we must proceed to the next order and estimate 



(R yk d k u y ) « (R yy ) [K(y)/y 2 ] 1/2 « c v (R yy ) [W(y)/y 2 ] 1/2 . (3.15) 
Using this this estimate we can proceed to conclude from Eq. (3.14) that 

±(P)(R yy )=g(R yy )W 1 / 2 ( y )/y, (3.16) 

where g is an empirical constant. Because equation (3.16) is crucial for the evaluation of 
the momentum flux in the buffer layer, we show in Fig. 3 (left panel) the quality of the 
closure approximation as obtained by numerical simulations. As one can see, the quality 
of the approximation is rather good within the buffer layer. 

This allows us to offer a final result for W(y) in which the polymer conformation tensor 
disappears altogether: 

W(y) ~ • ( 3 - 17 ) 

9 t 

This is an important result. On the one hand it shows how the Reynolds stress is reduced 
in the elastic layer, providing a direct mechanism to drag reduction. The Reynolds stress 
is the momentum transfer to the wall, and reducing it results in an increase of the mean 
momentum, in order to balance the constant rate of momentum production p'L. On the 
other hand, we see that this quantity is becoming of 0(y 2 ), justifying the neglect of 
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Figure 3. Left panehComparison between the term (Ryy){P)/T and (R yy )W 1 '" 2 /y in the nu- 
merical simulations. Right panehComparison of W as obtained by direct numerical simulations 
against the theoretical estimate y 2 {P) 2 /t 2 for two different values of r, namely r = 12.5 and 
t — 37.5. The estimate of g for the two cases gives g = 1.22 for r = 12.5 and g — 1.1 for 
r = 37.5. Thus, although the value of r changes by a factor 9, the combination of the Peterlin 
function with the t~ 2 term gives the same overall constant within 10 per cent error, showing 
the excellent quality of the prediction (3.17) 



Reynolds stresses in the simplified theory of L'vov et. al. (2003). In Fig. 3 (right panel) 
we compare the prediction of equation (3.17) against the direct numerical simulations. 
The quality of the results are extremely good, showing that the right physical behavior 
has been captured by employing the closure approximations. 

Neglecting then the first two terms in cq. (3.2), without the term p'y, and eq (3.11), 
and using the result (3.8) for the conformation tensor we end up with the estimate for 
the shear S(y) in the form 

Ag \/pT 



(3.18) 



v 2a v y 

Let us note that (3.18) is independent of r and it is consistent with the findings of 
L'vov et. al. (2003), i.e. it predicts a new logarithmic law of the wall which is the MDR 
asymptote. We note that the dissipative terms did not play any role in the viscoelastic 
derivation; the polymer terms took over, providing precisely the same role of the effective 
viscosity introduced in the simplified theory of L'vov et. al. (2003). 

In order to compute the value of A, we invoke the estimate of a in equation (3.11) as 
obtained in L'vov et. al. (2003), i.e. a = y^CN where y^ is the thickness of the viscous 
sub-layer, for the Newtonian fluid, expressed in the Prandtl units y + = y^p'L/u. Next, 
we compute the value of yF for which the MDR mean flow shear is equal to the viscous 
mean flow shear: 

Ag yfpJL 




which gives, in Prandtl units, 



Vf = -4= ■ ( 3 - 19 ) 



f 2a y 



For y = yF we expect that the energy dissipation due to viscous effect auK/yp is of the 
order of the energy dissipation due to polymers stretching, i.e. 
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Finally, by using the MDR estimate (3.18) and equations (3. 8), (3. 17, and (3.19), we can 
compute the value of A. It turns out that for the MDR slope Agj ^/2a y ~ y+(cAr/cy), 
in agreement with the previous findings of L'vov et. al. (2003). 



4. Saturation of drag reduction and cross over from the MDR to the 
Newtonian plug 

We expect a cross over from the MDR asymptote back to the Newtonian plug when 
the basic assumptions on the relative importance of the various terms in the balance 
equations lose their validity. Experimentally one sees that mean velocity profile follows 
the MDR up to some point yy after which it crosses over back to a logarithmic profile 
with the same slope as the Newtonian flow. To measure the quality of drag reduction 
one can introduce a dimensionless drag reduction parameter 

Q=^-l. (4.1) 

Vn 

The Newtonian flow is then a limiting case of the viscoclastic flow corresponding to 
Q — 0. There can be a number of reasons for the saturation of the drag reduction effect. 
One reason that can be important for large Reynolds numbers is that the concentration 
of the polymer is too small to provide the necessary increase in the effective viscosity as 
a function of the y + . This mechanism for the cross over back to the Newtonian plug was 
discussed in full detail in Benzi et. al. (2004). The main result of Benzi et. al. (2004) is 
that the drag reduction parameter Q is given by 

Q = a 3 c p Np CpSmall , Re large , (4-2) 

where c p is the number density (concentration in number of molecules per unit volume) 
of the polymer and N p is the degree of polymerization (the number of monomers per 
molecule). The parameter a is the linear scale of the monomer. This prediction was 
tested in Benzi et. al. (2004) by comparing with experiments with DNA as the drag 
reducing agent, with excellent agreement between theory and experiments. 

Here we address another mechanism for the saturation of drag reduction. We expect 
a cross over from the MDR asymptote back to the Newtonian plug when the basic 
assumptions on the relative importance of the various terms in the balance equations 
lose their validity, i.e. when (i) the turbulent momentum flux W becomes comparable to 
the total momentum flux p'L, or when (ii) the turbulent energy flux bK 3 ^ 2 jy becomes of 
the same order as the turbulent energy production WS. Using the estimates (3.17) and 
(3.18), one finds that both these conditions give the same cross-over point 



Note that f(y) = t/ (P(y)) is the effective non-linear polymer relaxation time. Therefore 
Eq. (4.3) can be also rewritten as 

S{y v )f{y v ) 1 • (4.4) 

In writing this equation we use the fact that the cross over point belongs also to the edge 
of the Newtonian plug where S(y) « ^p'L/y and that g <~ 1. The LHS of this equation 
is simply the local Weissenberg number (product of local mean shear and local effective 
polymer relaxation time). Thus, the cross-over to the Newtonian plug occurs at the point, 
where the local Weissenberg number decreases to ~ 1. We expect that this result remains 
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valid for any model of the elastic polymers, not only for the FENE-P model considered 
here. 

To understand how the cross-over point yv depends on the polymer concentration and 
other parameters, we need to estimate mean value of the Peterlin function (P) . Let us 
estimate the value of (P) as 

1 — a (it) 

where (R) = (R X x+Ryy + Rzz) ~ (R X x + 2R yy ) and a <~ (for simplicity we disregard 

p ). It follows from Eqs. (3.7)-(3.8), that 

(R xx ) ~ (Sf) 2 (R vv ) , 

and at the cross-over point (4.4) 

(Rxx) — (Ryy) , (R) — (Ryy) ■ 

The dependence of (R yy ) on y in the MDR region follows from Eqs. (3.11), (3.17), and 
(3.18): 

f 2 WS yVp^L 



(Ryy) : 

Then at the cross-over point y = yv'- 

(P>~ 



1 - ayv^/j/L/up 
Substituting this estimation into (4.3) gives the final result 

= CrWL 

U l + ap'Lr/vp V ; 

Here C is constant of the order of unity. Finally, introducing the dimensionless concen- 
tration of the polymers 

c P = , (4.6) 
a is 

we write the denominator in Eq. (4.5) as 

an' Lt 1 b'Lt , We 
1 + — = 1 + — = 1 + — ■ 

V p C p Vq Cp 

Here 

We ee tlL (4.7) 

is the (global) Weissenberg number. The final result for the dimensionless cross-over point 
jjy ee yv^/p'L/vQ assumes a very simple form: 

+ CWe 

»* " TTW5 ■ (48) 

This prediction can be put to direct test when c p is very large, or equivalently in the 
Oldroyd B model where P = 1 Bcnzi et. al. (2004). Indeed, in numerical simulations when 
the Weissenberg number was changed systematically, cf. Yu et. al. (2001), one observes 
the cross over to depend on We in a manner consistent with Eq. (4.8). The other limit 
when c p is very small is in agreement with the linear dependence on Cp predicted in Ecj. 
(4.2). 
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We can thus reach conclusions about the saturation of drag reduction in various limits 
of the experimental conditions, in agreement with experiments and simulations. 

5. Conclusions 

In this paper we showed how to understand some of the most common effects of dilute 
polymers solutions in wall bounded turbulent flows. Our starting point is the FENE-P 
model which has been shown to provide good qualitative and quantitative agreement 
with experimental findings. Within the context of the FENE-P model, we were able to 
derive in a controlled way the universal MDR profile that was first obtained theoretically 
in L'vov et. al. (2003). One of our main predictions is provided by Eq. (3.17) which tells 
us how the Reynold stress behaves as a function of the Peterlin function and the polymer 
relaxation time r. The predictions are in very good agreement with the results of DNS. 
In particular, Eq. (3.17) shows that the Reynolds stress is proportional to (_P) 2 t~ 2 , i.e. 
it becomes small for long polymers chains (r large) and large concentration. The velocity 
profile comes out independent on the concentration up to a point, yy, given by Eq. (4.8), 
where the von-Karman law is restored. In general our theoretical approach rationalizes 
and predicts a large number of qualitative and quantitative aspects of drag reduction as 
seen in laboratory experiments and in DNS of the FENE-P model. 

This work was supported in part by the European Commission under a TMR grant, 
the US-Israel Binational Science Foundation, and the Minerva Foundation, Munich, Ger- 
many. 
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